classdef Planet < handle
    %Planet Calculates orbit and angular momentum of planet
    %To solve this problem we consider Center of Mass frame.
    %The first approximation is ignoring the mass of planet as M(Sun)>>> m(Planet)
    
    properties
        a          %major axis in AU
        epsilon    %Elipticity
        dt         %timestep for RK4
        Period     %period of planet orbit in Year
        Tvec;	   %time vector
        GM=4*pi^2; %gravitational constant*Mass(Sun) in AU
        X0;
        X;
        Accel;	   %accelaration
    end
    
    methods

	% find the initial condition
        function calcX0(obj)
            x=obj.a*(1+obj.epsilon);
            v=sqrt(obj.GM*(2/x-1/obj.a));
	    obj.X0=[x;0;0;...
                0;v;0]; 
        end
        
        function calcOrbit(obj)
            obj.calcX0();
            P=ODESolver();
	    % A is -GM for Kepler's Orbit
            P.f=@(X)CentralPotential(X,-obj.GM);
            P.dt=obj.dt;
            P.T=obj.Period;
            P.X0=obj.X0;
            P.RK4();
            obj.X=P.X;
            obj.Tvec=P.Tvec;
	    % Calculate the accelaration
            r=sqrt( sum((P.X(1:3,:)).^2) );
            obj.Accel=[-obj.GM*P.X(1,:)./r.^3;...
                -obj.GM*P.X(2,:)./r.^3;...
                -obj.GM*P.X(3,:)./r.^3]
        end
        
        function plotOrbit(obj)
           figure

	   % plot Sun at origin
           plot3(0,0,0,'yo','linewidth',5)
           hold on

	   % Plot the planet 
           plot3(obj.X(1,:),obj.X(2,:),obj.X(3,:))
          
	   % Put accelaration vector at each point
           quiver3(obj.X(1,:),obj.X(2,:),obj.X(3,:),...
               obj.Accel(1,:),obj.Accel(2,:),obj.Accel(3,:))
           hold off
           axis equal
           xlabel('X (AU)')
           ylabel('Y (AU)')
           zlabel('Z (AU)')
        end
        
    end
    
end

